﻿#Make sure the paths are added first, such as by using
#sys.path.append(r'C:\My Code Projects\python-opencv\trunk\lib\python')

#The purpose of the demo will be to show a simple snake running in python.  Right now it is just showing the
#functionality of the snake module.

import Image
import ImageOps
import ImageFilter

import numpy as np
from pylab import *
from activeContour import snake as sn

x=np.zeros((201,201))
x[100,100]=1

y=sn.gaussianFilter(img=x,M=25,std=3)

figure
imshow(y)
title('Impulse response of the Gaussian filter');

x1=np.array([[1, 2,  3], [4, 5, 6]], np.float)
y1=sn.gradientMag(img=x1);
#expected answer
#array([[ 1.41421356,  2.23606798,  3.16227766],
#       [ 5.        ,  3.16227766,  3.16227766]])

#All of the vertex points
cp=np.array([33, 40, 64, 72, 70, 31, 27], np.float);
rp=np.array([32, 28, 29, 44, 71, 72, 53], np.float);



imR=numpy.asarray(ImageOps.grayscale(Image.open(r'C:\My Code Projects\activeContour\Jasons Code\circletest.png'))).astype(float);
y=sn.gaussianFilter(img=imR,M=15,std=3)
y2=sn.gradientMag(img=y);

#Take the difference of the points which will be used later 
rIdx=lambda x: x % p.shape[1];

#The difference between the points that define the snake contour will be used
#to compute the energy minimzation function
p=np.vstack([cp ,rp]);
pMagDerv=np.sqrt(np.dot(np.hstack([1, 1]),np.square(p-p[:,range(-1,p.shape[1]-1)])));
 
#compute the magnitude of the 2nd derivative 
tmpIdx=map(rIdx,range(1,p.shape[1]+1));
pMagDerv2=np.sqrt(np.dot(np.hstack([1, 1]),np.square(p[:,range(-1,p.shape[1]-1)]-2*p+p[:,tmpIdx])));

n=0
#for n in range(0,p.shape[1]):
print n

#we need to remove the indexes that correspond to the vertices which we will be testing new points
#the remaining idexes correspond to vertices that can be summed
partialSumIndexes=range(0,p.shape[1])
partialSumIndexes.remove(rIdx(n)); #n should never wrap, but might as well keep it since the code may get copied around
#we need to compute the partial sums without the current value of n
pMagDervPartialSum=np.sum(pMagDerv[partialSumIndexes])

partialSumIndexes.remove(rIdx(n-1));
partialSumIndexes.remove(rIdx(n+1));
pMagDerv2PartialSum=np.sum(pMagDerv2[partialSumIndexes])


#now we need to search a local region to see if we can minimize the energy function
#(newX,newY) = sn.minimizeE(rp,cp,n,gradMag,7,d);
neighborhoodSize=7     
#N = rp.size; #number of points in snake
halfSize = np.floor(neighborhoodSize/2);

#get local points around a snake vertex to be used for a local min search
(possR,possC) = mgrid[(rp[n]-halfSize):(rp[n]+halfSize+1),(cp[n]-halfSize):(cp[n]+halfSize+1)];

#we need to turn these into two column vectors that can be processed in parallel
possR=possR.flatten();
possC=possC.flatten();
pG=np.vstack([possC ,possR]);

Pk=np.vstack([possC ,possR]);
Pk_n1=tile(vstack(p[:,rIdx(n-1)]), (1, pG.shape[1]))  #n1 is neg 1
Pk_p1=tile(vstack(p[:,rIdx(n+1)]), (1, pG.shape[1]))  #p1 is positive 1

#all of the test cases for the average value

pGMagDerv=np.sqrt(np.dot(np.hstack([1, 1]),np.square(Pk-Pk_n1)));

pMagDervPartialSum=tile(pMagDervPartialSum, (1, pG.shape[1]))

d=(1.0/p.shape[1])*(pMagDervPartialSum+pGMagDerv);  #compute all of the Aaverage values


   
#we want the previous point, but the first point is special because it is duplicated 
#as the last point in the array so to get the previous point before the first we really need the true last 
#second to last point in the array.
if n==1:
    nAR = rp[N-1];
    nAC = cp[N-1];
else:
    nAR = rp[n-1];
    nAC = cp[n-1];

#get the next point
nBR = rp[n+1];
nBC = cp[n+1];

"""

imshow(y2)

hold('on');
plot(cp,rp,'r');
plot(cp,rp,'bx');
hold('off');
"""    


rp[n]=newX
cp[n]=newY
if n==0:  #we need to duplicate the end value
    rp[rp.size-1]=newX;
    cp[cp.size-1]=newY;
    
  



